#include "Substance.h"
#include <math.h>
#include <time.h>
#include <sys/timeb.h>
#include <stdio.h>

#include "common/common_defs.h"

///////Signum//////////
int sign(float a) {
    return (a > 0) - (a < 0);
}

float Particle::t;

int Component::Num_of_comp;
float Bound::t;
float Flow1::t_step;
IForce** Particle::ForceArr;


const float UGC = 8.314;

Component::Component(int molID, const MolInfo& mInfo, float Mass, int cID, float Temp) {
    Pos = mInfo.p;

    Vel = Vector3(rand_l(-100, 100), rand_l(-100, 100), rand_l(-100, 100));
    float absV = sqrtf(Temp * 3 * k_Boltz / Mass);
    Vel *= 1 / Vel.length() * rand_l(absV / 2, absV);

    CompID = cID;
    _molID = molID;

    M = Mass;

    avPeriod = 0;
    extForce = Vector3(0, 0, 0);

    Ek = Vel.length_2() * M / 2;
}
///////////////////////////////////////////////////////
void Component::InitComponent(const SubListItem& sub) {
    int NoC;
    int NoS;
    int NoC_1;
    Num_of_comp = sub.NUM_OF_COMPONENT;
    NoS = 1;

    ForceArr = new IForce *[Num_of_comp];
    for (NoC = 0; NoC < Num_of_comp; NoC++) {
        ForceArr[NoC] = new IForce[Num_of_comp];
    }
    for (NoC = 0; NoC < Num_of_comp; NoC++) {
        for (NoC_1 = 0; NoC_1 < Num_of_comp; NoC_1++) {
            ForceArr[NoC][NoC_1] = sub.ColFunc[NoC][NoC_1];
        }
    }
}

float Rand_rad(float R) {
    struct timeb tstruct;
    ftime(&tstruct);
    int time = tstruct.millitm;

    srand(time);
    float h = sin((float) rand());
    return h * R;
}

void Component::newPos() {
    Accel += extForce/M;

    avForce = dot(Accel, Vector3(1, 0, 0));

    Pos += Vel*t;
    Pos += Accel*(t * t);
    Vel += Accel*t;

    Ek = Vel.length_2() * M / 2;
    Accel = Vector3(0, 0, 0);
}
//////////////////////////////////////////
float Component::getAvForce() {
    return 0;
}

float Molecula::getAvForce() {
    float res = avForce;
    avForce = 0;

    return res;
}

void Component::normVel(float T, float T_f, int tau) {
    Vector3 v = Vel + Accel*t;
    float lVel = v.length();

    float c1 = 2 * M / 3.0 / k_Boltz;
    float T_fact = c1 * (lVel * lVel) / 2;

    if (T * 10 < T_fact) {
        Vel *= 1 / v.length() * sqrtf(T / c1);
        Accel *= 1 / v.length() * sqrtf(T / c1);
    } else {
        float alpha = sqrtf(1 - 1.0f / tau * (1 - T / (T_fact + 2)));
        Vel *= alpha;
    }

}

Molecula::Molecula(int idCh, const MolDscr& mDscr, const SubListItem& sub, float Temp) {
    ID = idCh;
    rank = mDscr.rank;
    compos = new Component*[rank];

    ///////////////////////////////////////////////
    Vector3 c_mass(0, 0, 0);
    Ek = 0;
    float Mass = 0;
    for (int i = 0; i < mDscr.rank; i++) {
        compos[i] = new Component(ID, mDscr.molSrc[i], sub.idName[mDscr.molSrc[i].cId].Mass, mDscr.molSrc[i].cId, Temp);
        c_mass += compos[i]->getPos();
        Mass += sub.idName[mDscr.molSrc[i].cId].Mass;
        Ek += compos[i]->Ek;
    }
    Ek = Ek / rank;
    printf(" Mmmm= %e \n", Ek);
    c_mass /= rank;
    Pos = c_mass;
    M = Mass;
}

TMolecula::TMolecula(int idCh, const TubeDscr * tDscr, const SubListItem& sub) {
    ID = idCh;
    rank = tDscr->rank;
    Temp = tDscr->Temp;
    rTau = tDscr->RelTime;
    OddEven = false;
    compos = new Component*[rank];
    distForce = new IForce[tDscr->NumOfPair];
    NofP = tDscr->NumOfPair;
    memcpy(distForce, tDscr->dCnt, NofP * sizeof(IForce));

    Vector3 c_mass(0, 0, 0);
    Ek = 0;
    float Mass = 0;
    for (int i = 0; i < tDscr->rank; i++) {
        compos[i] = new Component(ID, tDscr->molSrc[i], sub.idName[tDscr->molSrc[i].cId].Mass, tDscr->molSrc[i].cId, tDscr->Temp);
        c_mass += compos[i]->getPos();
        Mass += sub.idName[tDscr->molSrc[i].cId].Mass;
        Ek += compos[i]->Ek;
    }

    c_mass /= rank;
    Pos = c_mass;
    M = Mass;
}

void TMolecula::newPos() {
    float dist, vel, accel;
    const int RANK = this->NofP;

    for (int P = 0; P < RANK; P++) {
        int i = distForce[P]->I;
        int j = distForce[P]->J;

        Vector3 d = compos[i]->getPos() - compos[j]->getPos();
        Vector3 v = compos[i]->getVel() - compos[j]->getVel();

        dist = d.length();
        // vel = Vector3D::scalarProd(v.lenVector(), d.multVector(1 / dist));
        vel = dot(v, d/dist); // FIXME: either there _is_ or there _was_ a bug

        float M1 = compos[i]->getM();
        float M2 = compos[j]->getM();

        accel = Simple_Calc(distForce[P]->getForceInfo(), M1 * M2 / (M1 + M2), dist, t);

        if (accel > 1e+8)
            accel = 1e+5 / accel;

        Vector3 d_n = d/dist;
        Vector3 v_accel = d_n*accel;

        Vector3 v_accel1 = v_accel * (M2 / (M1 + M2));
        Vector3 v_accel2 = v_accel * (-M1 / (M1 + M2));

        compos[i]->addAccel(v_accel1);
        compos[j]->addAccel(v_accel2);
    }

    Pos = Vector3(0, 0, 0);
    Vel = Vector3(0, 0, 0);
    Ek = 0;
    float Tfact = 0;

    for (int i = 0; i < rank; i++) {
        Tfact += compos[i]->Ek * 2 / 3 / k_Boltz / rank * 2;
    }

    bool pr = false;
    for (int i = 21; i < rank - 21; i++) { // FIXME: WHY 21? T_T
        if (i % 2 == 0)
            pr = !pr;

        if (pr) {
            if (Tfact > 50)
                compos[i]->normVel(this->Temp, Tfact, this->rTau);

            compos[i]->newPos();
            Pos += compos[i]->getPos();
            Vel += compos[i]->getVel();
            Ek += compos[i]->Ek;
        } else
            compos[i]->Accel = Vector3(0, 0, 0);
    }
    Pos /= rank;
    Ek = Ek / rank * 2;

}

void Component::addExtForce(const Vector3& extF) {
    extForce += extF;
}

void Molecula::addExtForce(const Vector3& extF) {
    for (int i = 0; i < rank; i++)
        compos[i]->addExtForce(extF);

}

void Molecula::setPos(const Vector3& p, const Quaternion& o) {
    for (int i = 0; i < rank; i++) {
        Vector3 p1 = compos[i]->getPos() - Pos;
        Vector3 newP = twist(o, p1);
        compos[i]->setPos(p, Quaternion(0, newP));
        Ek = compos[i]->Ek;
    }
    Pos = p;

}
void TMolecula::setPos(const Vector3& p, const Quaternion& o) {
    for (int i = 0; i < rank; i++) {
        Vector3 p1 = compos[i]->getPos() - Pos;
        Vector3 newP = p1;
        newP = twist(o, newP) - p1;
        compos[i]->setPos(p, Quaternion(0, Vector3(0, 0, 0)));
        Ek = compos[i]->Ek;
    }
    Pos = p;
}

void Component::getPack(mpi_toSend& pack) {
    pack.sysID = _molID;
    pack.compId = CompID;
    pack.Pos = Pos;
    pack.Vel = Vel;
    pack.M = M;
    pack.Accel = Accel;
}

void Molecula::newPos() {
    avForce = 0;

    Pos = Vector3(0, 0, 0);
    Vel = Vector3(0, 0, 0);
    Ek = 0;
    for (int i = 0; i < rank; i++) {
        avForce += compos[i]->avForce * compos[i]->getM();
        compos[i]->newPos();
        Pos += compos[i]->getPos();
        Vel += compos[i]->getVel();
        Ek += compos[i]->Ek;
    }
    Ek = Ek / rank;
    Pos /= rank;
}

Flow1::Flow1(const SubListItem& s, MolDscr * mD, float P, float T, float R, float Rd)
: sub(s)
, pBuf(0)
, Temp(T)
, Ro(R * 1e-18)
, AvMass(0)
, MaxL(0)
, Rad(Rd)
{
    mDscr = new MolDscr[sub.nM];
    for (int i = 0; i < sub.nM; i++)
        mDscr[i] = mD[i];

    float tL = 0.35f;
    float tAvMass = 0;
    for (int i = 0; i < sub.nM; i++) {
        for (int j = 0; j < mDscr[i].rank; j++) {
            tL += 2 * mDscr[i].molSrc[j].p.length() / mDscr[i].rank;
            tAvMass += sub.idName[mDscr[i].molSrc[j].cId].Mass;
        }
        AvMass = AvMass + tAvMass * mDscr[i].Portion / 100;
        tAvMass = 0;
        if (tL > MaxL)
            MaxL = tL;
    }
    L = 5 * MaxL;
    L_TS = L;
    dRo = AvMass / (3.14 * Rad * Rad * L_TS * 2.56);
    MAX_LEN = 5 * floor(Ro / dRo);
    SP = 0;
    CUR = -1;
    NUM_P = 0;

    Prs = P * 3.14 * Rd * Rd / (Ro / dRo) * 1e-18;

    for (unsigned int i = 0; i < _NElem(streamMA); ++i)
        streamMA[i] = 0;

    maTop = 0;

    Pool = new Molecula*[MAX_LEN];
    for (int i = 0; i < MAX_LEN; i++)
        Pool[i] = 0;
}

void Flow1::normVel(int tau) {
    TempF = 0;
    float buf = 0;
    int ch = 0;

    for (int i = 0; i < CUR; i++)
        if (Pool[i] != 0) {
            buf += Pool[i]->Ek;
            ch++;
        }

    TempF = 2.0 / 3.0 * buf / ch / k_Boltz;

    for (int i = 0; i < CUR; i++)
        if (Pool[i] != 0)
            Pool[i]->normVel(Temp, TempF + 2, tau);
}

float Flow1::getTemp() {
    float result = 0;
    float buf = 0;
    int ch = 0;
    for (int i = 0; i < CUR; i++)
        if (Pool[i] != 0) {
            buf += Pool[i]->Ek;
            ch++;
        }

    result = 2.0 / 3.0 * buf / ch / k_Boltz;
    return result;
}

int Flow1::addPool(int idCh) {
    struct timeb tstruct;
    ftime(&tstruct);

    int time = 0;
    int time1 = 0;
    int time2 = 0;

    float c_chsr;
    int summ = 0;

    float factRo = NUM_P * dRo;
    Diver = floor((Ro - factRo) / dRo);

    delete[] pBuf;
    pBuf = new Molecula*[Diver];

    for (int i = 0; i < Diver;) {
        ftime(&tstruct);
        time += tstruct.millitm + 1;
        c_chsr = 50 * sin((float) time) + 50;
        summ = 0;

        for (int j = 0; j < sub.nM; j++) {
            summ += floor(mDscr[j].Portion);
            if (summ / (int) ceil(c_chsr) != 0) { // FIXME: why we divide when checking "!= 0"?
                while (Pool[SP] != 0)
                    SP++;

                Pool[SP] = new Molecula(idCh, mDscr[j], sub, Temp);

                ++idCh;

                pBuf[i] = Pool[SP];

                //////////////////////////////////
                ftime(&tstruct);
                time += tstruct.millitm;
                ftime(&tstruct);
                time1 += tstruct.millitm;
                ftime(&tstruct);
                time2 += tstruct.millitm;

                // FIXME: why the hell we use such crazy distribution?
//                Vector3 o = Vector3(sin((float) time) / 1000.0, sin((float) time1) / 1000.0, sin((float) time2) / 1000.0);
//                float tw = sin((float) time1 * (float) time / 1000.0);
//                Quaternion t1 = Quaternion(tw, o);
//                t1.setNorm(1.0f);
                Vector3 o(sin((float) time), sin((float) time1), sin((float) time2));
                float tw = sqrtf(1 - o.length_2()); // yes i am too lazy to implement norm()/setNorm() just for using it here just once
                Quaternion t1 = Quaternion(tw, o);
                ftime(&tstruct);
                time += tstruct.millitm;
                float y = Rad * sin((float) time / 1000.0);
                ftime(&tstruct);
                time1 += tstruct.millitm;
                float z = Rad * cos((float) time1 + time / 1000.0);
                ftime(&tstruct);
                time += tstruct.millitm;
                Vector3 p((L_TS / 6 * sin((float) time) + L_TS / 6), y, z);

                Pool[SP]->setPos(p, t1);
                Pool[SP]->addExtForce(Vector3(Prs, 0, 0));

                ///////////////////////////////////

                if (SP > CUR)
                    CUR++;

                NUM_P++;

                i++;
                break;
            }
        }
    }

    float Ek_n = Temp * 3 / 2 * 1.38e-14;

    float ax_pos;
    int dst_ch = 0;
    for (int i = 0; i < CUR; i++) {
        if (Pool[i] != 0) {
            ax_pos = Pool[i]->getPos().x;
            if (ax_pos > (L_TS + 0.4) && Pool[i]->Ek < Ek_n * 3) {
                if (i == CUR)
                    CUR--;
                if (i < SP)
                    SP = i;
                NUM_P--;
                dst_ch++;
                Pool[i] = 0;        // FIXME: memleak?
            }
        }
    }
    streamMA[maTop % _NElem(streamMA)] = dst_ch;
    maTop++;

    return Diver;
}

float Flow1::getStream() {
    float summ = 0;
    for (unsigned int i = 0; i < _NElem(streamMA); i++) {
        summ += streamMA[i];
    }
    summ /= (float)_NElem(streamMA);

    summ = summ / t_step * AvMass;

    return summ;
}

Flow1::~Flow1() {
    delete[] Pool;      // FIXME: how about deleting array contents?
    delete[] pBuf;
    delete mDscr;
}

RBound::RBound(const Vector3& pos, const Vector3& axis, float rad, float h) {
    Axis = axis;
    Pos = pos;
    param1 = rad;
    param2 = h;
}

bool RBound::belon(const Vector3& lV, const Vector3& rV) {
    Vector3 coordL = lV - Axis*dot(lV, Axis);
    Vector3 coordR = rV - Axis*dot(rV, Axis);

    float coord[] = { coordL.y, coordL.z, coordR.y, coordR.z };
    Vector3 vTop(0, coord[0], coord[0]);
    float rTop = vTop.length();
    bool pr;

    if (rTop <= param1)
        pr = true;
    else
        pr = false;

    for (int i = 0; i < 4; i++)
        for (int j = 0; j < 4; j++) {
            if (i != 0 && j != 0)
                vTop = Vector3(0, coord[i], coord[j]);
            rTop = vTop.length();       // FIXME: should be under previous if block?
            if (rTop >= param1 && pr)
                return true;
            if (rTop <= param1 && (!pr))
                return true;
        }

    return 0;
}

PlBound::PlBound(const Vector3& pos, const Vector3& norm, float RadInt, float RadExt, int apc) {
    Pos = pos;
    Axis = norm;
    param1 = RadInt;
    param2 = RadExt;
    AvPrsCh = apc;

    prsBuf = new float[AvPrsCh];
    for (int i = 0; i < AvPrsCh; i++)
        prsBuf[i] = 0;

    prs = 0;

    Sqr = 3.1415 * RadExt * RadExt - 3.1415 * RadInt * RadInt;
    bufP = 0;
}

bool PlBound::belon(const Vector3& lV, const Vector3& rV) {
    bool pr = false;

    Vector3 n1 = Pos - lV;
    Vector3 n2 = Pos - rV;
    float d1 = dot(n1, Axis);
    float d2 = dot(n2, Axis);

    Vector3 tan1 = cross(n1, Axis);
    Vector3 tan2 = cross(n2, Axis);

    if ((tan1.length() > param1 && tan1.length() < param2) || (tan2.length() > param1 && tan2.length() < param2)) {
        pr = true;
    }

    if (sign(d1) == sign(d2)) {
        pr = false;
    }

    return pr;
}

void PlBound::bCtrl(Component* p) {
    IFr a(0, 0, 0.2, 1e-12, 'L');
    IFr force(0, 0, 0.1, 1e-12, 'L');

    Vector3 p_c = p->getPos() - Pos;
    float dist = fabsf(dot(Axis, p_c));

    float R = cross(p_c, Axis).length();

    bool pr = false;
    if (R > param1 && R < param2) {
        pr = true;
    }

    if (dist < a.Potential(0, 'R') * 1.01f && pr) {
        float accel = Simple_Calc(a.getForceInfo(), p->getM(), dist, t);
        Vector3 v_accel = Axis*(dot(p_c, Axis) * accel);
        p->addAccel(v_accel);

        prs += dot(Axis, v_accel) * p->getM() / Sqr;
    }
    Vector3 np_c = p_c - Axis*dist;
    np_c = np_c*(1 / np_c.length() * param1);
    p_c -= np_c;

    if (param1 != 0) {
        if (dist < a.Potential(0, 'R') * 1.01) {
            float rDist = p_c.length();
            if (rDist < a.Potential(0, 'R') * 1.01) {
                float accel = Simple_Calc(force.getForceInfo(), p->getM(), rDist, t);
                Vector3 v_accel = p_c*(1 / p_c.length() * accel);
                p->addAccel(v_accel);
            }
        }
    }
}

void PlBound::CalcPrs() {
    prsBuf[bufP % AvPrsCh] = prs;
    bufP++;
    prs = 0;
}

float PlBound::getPrs() {
    float res = 0;
    for (int i = 0; i < AvPrsCh; i++)
        res += prsBuf[i];

    return res / AvPrsCh * 10e+18;
}

void RBound::bCtrl(Component* p) {
    IFr a(0, 0, 0.1, 1e-12, 'L');

    Vector3 p_c = p->getPos();
    float scal = dot(p_c, Axis);
    Vector3 aV1 = Axis*scal;
    p_c -= aV1;
    float len = p_c.length();
    float dist = param1 - len;

    if (dist < a.Potential(0, 'R') * 1.0f) {
        Vector3 cross_n = p_c/len;
        float accel = Simple_Calc(a.getForceInfo(), p->getM(), dist, t);
        Vector3 v_accel = -accel*cross_n;

        p->addAccel(v_accel);
    }
}

void Flow1::prsControl(float g) {

    float diver = -10e-13 - g;
    diver = diver * 3.14 * Rad * Rad * 1.96;

    float d_force = diver / CUR;

    for (int i = 0; i < CUR; i++) {
        if (Pool[i] != 0)
            Pool[i]->addExtForce(Vector3(-d_force, 0, 0));
    }
}

float Flow1::getPrs() {
    float avF = 0;
    int ch = 0;
    for (int i = 0; i < CUR; i++) {
        if (Pool[i] != 0)
            if (Pool[i]->getPos().x > L_TS / 2) {
                avF += Pool[i]->getAvForce();
                ch++;
            }
    }

    avF = avF / ch / 3.14 / Rad / Rad / 1.96;

    return avF;
}

void Component::setPos(const Vector3& p1, const Quaternion& q) {
    Pos += p1 + q.toVector();
}

void Molecula::normVel(float T, float Tf, int tau) {
    for (int i = 0; i < rank; i++)
        compos[i]->normVel(T, Tf, tau);
}
